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ABSTRACT 

Lambda Cold Dark Matter (ACDM) is now the standard theory of structure formation in the 
Universe. We present the first results from the new Bolshoi dissipationless cosmological ACDM 
simulation that uses cosmological parameters favored by current observations. The Bolshoi sim- 
ulation was run in a volume 250 h~ l Mpc on a side using ~8 billion particles with mass and 
force resolution adequate to follow subhalos down to the completeness limit of V clrc — 50 km s -1 
maximum circular velocity. Using merger trees derived from analysis of 180 stored time-steps 
we find the circular velocities of satellites before they fall into their host halos. Using excellent 
statistics of halos and subhalos (~10 million at every moment and ^50 million over the whole 
history) we present accurate approximations for statistics such as the halo mass function, the 
concentrations for distinct halos and subhalos, abundance of halos as a function of their circular 
velocity, the abundance and the spatial distribution of subhalos. We find that at high redshifts 
the concentration falls to a minimum value of about 4.0 and then rises for higher values of halo 
mass, a new result. We present approximations for the velocity and mass functions of distinct 
halos as a function of redshift. We find that while the Sheth-Tormen approximation for the 
mass function of halos found by spherical overdensity is quite accurate at low redshifts, the ST 
formula over-predicts the abundance of halos by nearly an order of magnitude by z = 10. We find 
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that the number of subhalos scales with the circular velocity of the host halo as U host , and that 
subhalos have nearly the same radial distribution as dark matter particles at radii 0.3-2 times 
the host halo virial radius. The subhalo velocity function N(> V su b) scales as V^.. Combining 
the results of Bolshoi and Via Lactea-II simulations, we find that inside the virial radius of halos 
with V^i rc = 200 km s _1 the number of satellites is N(> U su b) = (U su b/58 km s -1 )~ 3 for satellite 
circular velocities in the range 4 km s _1 < V^ u b < 150 km s . 



Subject headings: cosmology: theory - 
numerical 
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1. Introduction 

The Lambda Cold Dark Matter (ACDM) model 
is the standard modern theoretical framework for 
understa nding the formation of structure in the 
universe ( Dunklev et al.l [20091) . With initial con- 
ditions consisting of a nearly scale-free spectrum 
of Gaussian fluctuations as predicted by cosmic 
inflation, and with cosmological parameters de- 
termined from observations, ACDM makes de- 
tailed predictions for the hierarchical gravita- 



tional growth of structure. For the past sev- 
eral years, the best large simulation for com- 
parison with galaxy surveys has been th e Mil- 



lennium Simulation (Springc l et al. 20051 MS 



I). Here we present the first results from a new 
large cosmological simulation, which we are call- 
ing the Bolshoi simulation ("Bolshoi" is the Rus- 
sian word for "big'' Q ). Bolsho i has nearly an 



1 "Bolshoi" can be translated as (1) big or large (2) great (3) 
important (4) grown-up. The Bolshoi Ballet performs in 
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order of magnitude better mass and force reso- 
lution than the Millennium Run. The Millen- 
nium Run used the first-year (WMAPf) cosmo- 
logical parameters from the Wilkinson Microwav e 
Anisotropy Probe satellite (jSpergel et al, 2003). 
These parameters are now known to be incon- 
sistent with modern measurements of the cos- 
mological parameters. Th e Bolshoi simulation 
used the lates t WM AP5 (iHinshaw et al.l [2009t 



Komatsu et all 120091: fPunklev et al.l l2009t ) and 



dius and r s is the characteristic radius where the 
log-log slope of the density is equal to —2. Details 
are given in §3 and in Appendix B. 

One of the main goals of this paper is to pro- 
vide the basic statistics of halos selected by the 
maximum circular velocity ( V^i rc ) of each halo. 
There are advantages to using the maximum cir- 
cular velocity as compared with the virial mass. 
The virial mass is a well defined quantity for dis- 
tinct halos (those that are not subhalos), but it 
is ambiguous for subhalos. It strongly depends 
on how a particular halofinder defines the trun- 
cation radius and removes unbound particles. It 



1984t iBlumenthal et al.lll984f) soon led to the first 



CDM A^-body cosmologi cal simulations ([Melott et al 



1983c iDavis et al.lll985j ). Ever since then, such 



WMAP7 (jJarosik et al.l 120101 ) parameters, which 
are also consistent with other recent observational 
data. 

The invention of Cold Dark Matter ( Primack fc Blumenfehadi depends on the distance to the center of the 

host halo because of tidal stripping. Instead, the 
circular velocity is less prone to those complica- 
tions. Even for distinct halos the virial mass is an 
inconvenient property because there are different 
definitions of "virial mass" . None of them is better 
than the other and different research groups prefer 
to use their own definition. This causes confusion 
in the community and makes comparison of re- 
sults less accurate. The main motivation for using 
Vdrc is that it is more closely related to the prop- 
erties of the central regions of halos and, thus, to 
galaxies hosted by those halos. For example, for a 
Milky- Way type halo the radius of the maximum 
circular velocity is about 40 kpc (and the circular 
velocity is nearly the same at 20 kpc), while the 
virial radius is about 300 kpc. As an indication 
that circular velocities are a better quantity for 
describing halos, we find that most statistics look 
very simple when we use circular velocities: they 
arc cither pure power-laws (abundance of subha- 
los inside distinct halos) or power-laws with nearly 
exponential cutoffs (abundance of distinct halos). 

This paper is organized as follows. In §2 we give 
the essential features of the Bolshoi simulation. 
Section 3 describes the halo identification algo- 
rithm used in our analysis. In §4 we present results 
on masses and concentrations of distinct halos. 
Here we also present relations between V c i rc and 
M v i r . The halo velocity function is presented in §5. 
Estimates of the Tully-Fish er relation are given in 
Truiillo-Gomez et al.l < 201o[ l . where we present de- 



simulations have been essential in order to calcu- 
late the predictions of CDM on scales where struc- 
ture has formed, since the nonlinear processes of 
structure formation cannot be fully described by 
analytic calculations. For example, one of the 
first large simulatio ns of the ACDM cosmology 
dKlvpin et al.lll996h showed that the dark matter 
autocorrelation function is much larger than the 
observed galaxy autocorrelation function on scales 
of ~ 1 Mpc, so "scale-dependent anti-biasing" 
was required for ACDM to match the observed 
distribution of galaxies. Subsequent simulations 
with resolution adequate to id entify the dark mat- 
ter halos that ho st galaxies ( Jenkins et al. 1998; 
Colin et al.l [T999) demonstrated that the required 



destruction of dark matter halos in dense regions 
does indeed occur. 

./V-body simulations have been essential for 
determining the properties of dark matter ha- 
los. It turned out that dark matter halos of 
all masses typically have a similar radial pro- 
file, wh ich can be approxi mated by the NFW 
profile ( Navarro et al.l Il996l ). Simulations were 
also crucial for determining the dependence of 
halo concentration c v - n - = R v ^/r s and halo shape 
on halo mass and redshift ( Bullock et al.l 2001 



Zhao- et al.l 120031: lAlkood et all 120061: iNeto etal 
20071 iMaccio et al.1 120081 ) , and also for determin- 



ing the dependence of the concen tration of halos 



on their mass accreti on history ([Wechsler et al 



2002: IZhao et al.1 12009;). Here i? v i r is the virial ra- 



Bolshoi Theater in Moscow. 



tailed discussions of numerous issues related with 
the procedure of assigning luminosities to halos in 
simulations and confront results with the observed 
galaxy distribution. In §6 we give statistics of the 
abundance of subhalos. The number density pro- 
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Table 1 
cosmological parameters 



Hubble h 



tilt n 



ref. 



WMAP5 

WMAP5+BAO+SN 
X-ray clusters 

X-ray clusters+WMAP5 
X-ray clusters+WMAP5 

+SN+BAO 
maxBCG + WMAP5 
WMAP7+BAO+H 
Bolshoi simulation 
Millennium simulations 
Via Lactea-II simulation 



0.719 i 



1+0.026 
-0.027 

0.701 ±0.013 
0.715 ±0.012 



(0.719) 
(0.72 ±0.08) 



(0.70) 
0.704 
0.70 
0.73 
0.73 



+0.013 



0.258 ± 0.030 
0.279 ±0.013 
0.260 ±0.012 



-,+0.03 
-0.02 

0.269 ±0.016 



0.30^ 



0.265 ±0.016 
0.273 ±0.014 
0.270 
0.250 
0.238 



+0.014 
0.015 
0.1114 



0.963 
0.960_^ 
(0.95) " 

(0.963) 
(0.95) 



(0.96) 

0.963 ±0.012 

0.95 

1.00 

0.951 



0.796 ±0.0326 

0.817 ±0.026 

0.786 ± 0.011 

±0.020 a 

S5+ 04 
u -°°-0.02 

0.82 ±0.03 

0.807 ±0.020 

0.809 ± 0.024 

0.82 

0.90 

0.74 



Dunklcy ct al. (2009) 
Dunklev et al. (2009) 
Vikhlinin et al. (2009) 



Henry et al. (2009) 
Mantz et al. (2008) 



Rozo et al. (2009) 
Jarosik et al. (2010) 



Springcl ct al. (2005) 
Diemand et al. (2008) 



Note. — Values in parentheses are priors 
a Systematic error 



files of subhalos are presented in §7. Section §8 
gives a short summary of our results. Appendix 
A describes details of the halo identification pro- 
cedure. Appendix B collects useful approximation 
formulas. Appendix C compares masses of halos 
found with two different halo- finders: the friends- 
of-friends algorithm and the spherical overdensity 
halo-finder used in this paper. 

2. Cosmological parameters and Simula- 
tions 

The Bolshoi simulation was run with the 
cosmological parameters listed in Table [TJ to- 
gether with fibar = 0.0469, n = 0.95. As 
Table Q] shows, these parameters are compati- 
ble with the WMA P seven-year data (WMAP7) 
(jjarosik et al.ll2010h results, with the WMAP five- 
year data (WMAP5), with WMAP 5 combined 
with Baryon Acoustic Osci l lations and Supernova 
data dHinshaw et all [20091 iKomatsu etall l2009t 



Dunklev et all l2009h . The parameters used for 



the Bolshoi simulation are also compatible with 
the recent constraints fro m the Chandra X-ray 
cluster cosmology project ( Vikhlinin et al.l 12009) 



and other recent X-ray cluster studies. The Bol- 
shoi parameters are in excellent agreement with 
the SPSS maxBCG+WMA P5 cosmological pa- 
rameters ( Rozo et al. 20091) . The optical cluster 



measurements of the SDSS maxBCG cluster cat- 
alog are fully consistent with the WMAP5 data, 
and the joint maxBCG+WMAP5 analysis quoted 
in Table [1] reduces the errors. Figure [T] shows 
current observational constraints on the Dm and 
ct 8 parameters. It shows graphically that Bolshoi 
agrees with the recent constraints while Millen- 
nium is far outside them. 

The Millennium simulation ([Springel et al 



abundance and weak gravitational lensing mass 



2005, MS-I) has been the basis for many stud- 
ies of the distribution and statistical properties 
of dark matter halos and for semi-analytic mod- 
els of the evolving galaxy population. However, 
it is important to appreciate that this simula- 
tion a nd the more recent Millen nium-II simu- 
lation (jBovlan-Kolchin et al.l l2009l MS-II) used 
the first-year (WMAP1) cosmological parameters, 
which are rather different from the current param- 
eters summarized in Table [TJ The main difference 
is that the Millennium simulations used a sub- 
stantially larger amplitude of perturbations than 
Bolshoi. Formally, the value of as used in the Mil- 
lennium simulations is more than 3<r away from 
the WMAP5+BAO±SN value and nearly 4cr away 
from the WMAP7+BAO+H value. However, the 
difference is even larger on galaxy scales because 
the Millennium simulations also used a larger tilt 
n = 1 for the power spectrum. Figure [2] shows the 
linear power spectra of the Bolshoi and Millen- 
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simulations and testing particular characteristics. 



i Millennium 



Mantz et al. 2008 

Henry et al. 2009 

Vikhlinin et al. 2009 

Rozo et al. 2010 




olshoi \ 



0.92 
0.90 - 
88 
0.86 
0.84 
0.82 
0.80 - 
0.78 
0.76 
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Fig. 1. — Optical and X-ray cluster abundance 
plus WMAP constraints on as and f2jw. Con- 
tours show 68% confidence regions for a joint 
WMAP 5 and cluster abundance analysis assum- 
ing a flat ACDM cosmology. The shaded re- 
gion is the SDSS optical max BCG cluster abun- 
dance + WMAP5 analysis from lRozo et all (|2009h 
(which is also the source of this figure). The 
X-ray + WMAP 5 constraints are from several 
sourc es: the low-redshif t cluster luminosity func- 
tion ( Mantz et al. 20081). th e cluster temperature 
function (Hen ry et al. 112009ft. and the cluster mass 



1.8 - 



function ( Vikhlinin et al. 1 12009). All four recent 



studies are in excellent agreement with each other 
and with the Bolshoi cosmological parameters. 
The Millennium I and II cosmological parameters 
are far outside these constraints. 



nium simulations. Because of the large difference 
in the amplitude, it is not surprising that the Mil- 
lennium simulations over-predict the abundance 
of galax y-size halos. The Sheth -Tormen approx- 



imation (jSheth fc Tormenl 120021 ) gives a factor of 
1.3-1.7 more Mvir ~ lO 12 /!r 1 M halos at z = 2-3 
for Millennium as co mpared with Bolshoi, w hich 
is a large difference. lAngulo fc White (2009) ar- 
gue that cosmological N-body simulations can be 
re-scaled by certain approximations or by other 
means. However, the accuracy of those re-scalings 
cannot be estimated without running accurate 
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Fig. 2. — Bottom: Linear power spectra of the Bol- 
shoi and Millennium simulations at redshift zero. 
Top: Ratio of the spectra. The Millennium simu- 
lations have substantially larger amplitude of per- 
turbations on all scales, resulting in overpredic- 
tion of the number of galaxy-size halos at high 
redshifts. 



Table 2: Parameters of Bolshoi simulation 



Parameter 



Value 



Box size (h~ L Mpc) 


250 


Number of particles 


2048 3 


Mass resolution (/("'Mq) 


1.35 x 10 8 


Force resolution 




(/i _1 kpc, physical) 


1.0 


Initial redshift 


80 


Zero-level mesh 


256 3 


Maximum number of 




refinement levels 


10 


Zero-level time-step Aa 


(2-3) x 10~ 3 


Maximum number of time-steps 


~ 400, 000 


Maximum displacement 


0.10 


per time-step 


(cell units) 



The Bolshoi simulation uses a computational 
box 250/i _1 Mpc across and 2048 3 rj 8 billion par- 
ticles, which gives a mass resolution (one parti- 
cle mass) of mi = 1.35 x 10 8 h' 1 M G . The force 
resolution (smallest cell size) is physical (proper) 
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1 h~ kpc (see below for details). For compari- 
son, the Millennium-I simulation had force reso- 
lution (Plummer softening length) 5 /i _1 kpc and 
the Millennium-II simulation had 1 ft, kpc. Ta- 
ble [5] gives a short summary of various numerical 
parameters of the Bolshoi simulation. 

The Bolshoi simulation was performed with the 
Adaptivc-Rcfinement-Tree (ART) code, which is 
an Adaptive-Mesh-Refinement (AMR) type code. 
A detailed descriptio n of the code is giv en in 
Kravtsov et"aH (|l997t ) and iKravtsovl (|l999t) . The 



code was parallelize d using MPI libraries an d 
OpenMP directives (|Gottloeber fc Klvpinl [iooih . 



Details of the time-stepping algorithm and com- 
parison with GADGET and PKDGRAV codes are 
given in lKlypin et al. (2009). Here we give a short 



outline of the code and present details specific for 
Bolshoi. 

The ART code starts with a homogeneous mesh 
covering the whole cubic computational domain. 
For Bolshoi we use a 256 3 mesh. The Cloud-In- 
Cell (CIC) method is used to obtain the density on 
the mesh. The Poisson equation is solved on the 
mesh with the FFT method with periodic bound- 
ary conditions. The ART code increases the force 
resolution by splitting individual cubic cells into 
2x2x2 cells with each new cell having half the 
size of its parent. This is done for every cell if the 
density of the cell exceeds some specified thresh- 
old. The value of the threshold varies with the 
level of refinement and with the redshift. Once 
the hierarchy of refinement cells is constructed, the 
Poisson equation is solved on each refinement level 
using the Successive OverRelaxation (SOR) tech- 
nique with red-black alternations. Boundary con- 
ditions arc taken from the one-level coarser grid. 
The initial guess for the gravitational potential is 
taken from the previous time-step whenever pos- 
sible. 

We use the on-line Code for Anisotr opies in 



the M icrowave Background of lLewis et al.l (CAMB 



2000 jl to generate the power spectrum of cosmo- 



logical perturbation s. The code is based on th e 
CMBFAST code bv lSeliak fc Zaldarriaeal (|l996h . 



At early moments of evolution, when the am- 
plitude of perturbations was small, the refinement 
thresholds were chosen in a way that allows unim- 
peded growth of even the shortest perturbations, 



http: //lambda. gsfc .nasa .gov/toolbox/tb_camb_f orm. cfm 








k (Mpc/h) 

Fig. 3. — Growth of the power spectrum of 
perturbations at early stages of evolution. The 
full curves show A 2 = fc 3 P(fc)/27r 2 at redshifts 
z = 11,20,53 (from top to bottom). The dashed 
curves show the linear power spectrum. The ver- 
tical line is the Nyquist frequency of particles. 

with wavelengths close to the Nyquist frequency. 
Ideally, the distance between particles should be 
at least two cell sizes: at that separation the force 
of gravity is Newtonian. In setting parameters for 
Bolshoi we came close to this condition: we used 
a density threshold of 0.6 particles per cell, which 
resulted in effectively resolving the whole compu- 
tational volume down to 4 levels of refinement or, 
equivalently, to a 4096 3 mesh. This condition was 
kept until redshift z = 11. As the perturbations 
grow, get nonlinear, and collapse, it becomes pro- 
hibitively expensive (memory consuming) and not 
necessary to keep this strict refinement condition. 
We gradually increase the threshold for the 4-th 
refinement level: 0.8 particles till z = 9, 1 parti- 
cle till z = 7, 3 particles till z = 1.5, and there- 
after we had 5 particles. The same increase of the 
thresholds was done for higher refinement levels, 
for which we started with 2 particles at high z and 
ended with 5 particles at z = 0. 

Figure [3] shows the power spectrum of pertur- 
bations at redshifts z = 11, 20, 53 and compares it 
with the linear theory. A 4096 3 mesh was used to 



5 



estimate P(k) in the simulation, which may un- 
derestimate the real spectrum by 3-5 percent at 
the Nyquist frequency due to the finite smoothing 
of the density field produced by the Cloud-In-Cell 
density assignment. Results indicate that the code 
was evolving the system as it was expected to: lin- 
ear growth of small perturbations at early stages 
and at large scales with no suppression at high 
frequencies. 

The number of refinement levels, and thus the 
force resolution, change with time. For example, 
there were 8 levels at z = 10, and 9 at z = 5, 
which gives the proper resolution of 0.4 h~ 1 kpc. 
At z = 1 the tenth level was open with the proper 
resolution of 0.5 /i _1 kpc. However, when a level of 
refinement opens, it contains only a small fraction 
of volume and particles. Only somewhat later does 
the number of particles on that level become sub- 
stantial. This is why the evolution of the force res- 
olution is consistent with nearly constant proper 
force resolution of 1 ft, _1 kpc from z — 20 to z = 0. 

The ART code uses the expansion parameter 
a = (1 + as the time variable. Each parti- 
cle moves with its own time-step which depends 
on the refinement level. The zero-level defines the 
maximum value Aa of the stepping of the expan- 
sion parameter. For Bolshoi Aa = 2 x 10~ 3 for 
a < 0.8 and Aa = 3 x 10~ 3 for later moments. 
The time-step decreases by a factor of two from 
one level of refinement to the next. There were 
ten levels of refinement at the later stages of evo- 
lution (z < 1). This gives the effective number of 
400,000 time-steps. 

The initial conditions for the simulation were 
created using the Zeldovich approxi mation with 



tion. For this particle the density perturbation 
is 6.5 x 0.0826 = 0.55; still below 1.0. The most 
dense 100 particles are expected to have a density 
contrast 5.76 x 0.0856 = 0.49, low enough for the 
Zeldovich approximation still to be accurate. 

The Bolshoi simulation was run at the NASA 
Ames Research Center on the Pleiades supercom- 
puter. It used 13,824 cores (1728 MPI tasks 
each having 8 OpenMP threads) and a cumulative 
13 Tb of RAM. We saved 180 snapshots for sub- 
sequent analysis. The total number of files saved 
in different formats is about 600,000, which use 
100 Tb of disk space. 

For some comparisons we also use a cata- 
log of halos fo r the Via Lactea-II simulation 



( Diemand et al. 20081 ). This simulation was 



run for one isolated halo with maximum cir- 
cular velocity V C i IC = 201 km s -1 . Using ap- 
proximatio ns for the dar k matt er profile pro- 
vided by iDiemand et~aH (|2008l) . we estimate 
the virial mass and radius of the halo to be 
M vir = 1.3 x 10 12 h" 1 M Q and R vil = 226/T 1 kpc. 
Here we use the top-hat model with cosmologi- 
cal constant A to estimate the virial radius, as 
explained in the next section. The Via Lactea-II 
simulation uses slightly different cosmological pa- 
rameters as compared with Bolshoi (see Table 1). 
In particular, the amplitude of perturbations is 10 
percent lower: <rg = 0.74. 

3. Halo identification 

Let us start with some definitions. A distinct 
halo is a halo that does not "belong" to another 
halo: its center is not inside of a sphere with a 



. , , , . , nr, — . — r~cn j — radius equal to the virial radius of a larger halo. 

particles placed m a homogeneous grid (Klypm & bhandarin , . ..... . . . 

i nool It^i — o tt ij. — IN nrnK ni^ subhalo is a halo whose center is msidc tfic virial 

1983; Klypm &: Holtzmanlll997l l. Bolshoi starts at 



^init = 80 when the rms density fluctuation Ap/p 
in the computational box with 2048 3 particles was 
equal to Ap/p = 0.0826. (The Nyquist frequency 
defines the upper cutoff of the spectrum, with the 
low cutoff being the fundamental mode.) The ini- 
tialization code uses the Zeldovich approximation, 
which provides accurate results only if the den- 
sity perturbation is less than unity. This should 
be valid at every point in the volume. In prac- 
tice, the fluctuations must be even smaller than 
unity for high accuracy. With 2048 3 independent 
realizations of the density, one expects to find 
one particle in the box to have a 6.5<r fluctua- 



radius of a larger distinct halo. Note that distinct 
halos may overlap: the same particle may belong 
to (be inside the virial radius of) more than one 
halo. We call a halo isolated if there is no larger 
halo within twice its virial radius. In some cases 
we study very isolated halos with no larger halo 
within three times its virial radius. 

We use the Bound-Density-Maxi ma (BDM) al- 

gorit hm to identify halos in Bolshoi ([Klypin fc Holtzman 
19971 ). App endix A gives s ome o f the details of the 



halofinder. iKnebe et al. ( 2011 ) present detailed 
comparison of BDM with other halofinders and 
show results of different tests. The code locates 
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maxima of density in the distribution of particles, 
removes unbound particles, and provides several 
statistics for halos including virial mass and ra- 
dius, and maximum circular velocity 



GM{< r) 



(1) 



Throughout this paper we will use V^h- C and the 
term circular velocity to mean maximum circular 
velocity over all radii r. When a halo evolves over 
time, its V c ; rc may also evolve. This is especially 
important for subhalos, which can be significantly 
tidally stripped and can reduce their V^rc over 
time. Thus, we distinguish instantaneous V^i rc and 
the peak circular velocity over halo's history. 



We use the virial mass definition M v 



that 



follows from the top-hat model in an expanding 
Universe with a cosmological constant. We de- 
fine the virial radius i? v ir of halos as the radius 
within which the mean density is the virial over- 
density times the mean universal matter density 
p m = f^M/fcrit at that redshift. Thus, the virial 
mass is given by 



An 3 

M v i r = -^-AvirPmi? v i; 



(2) 



For our set of cosmological parameters, at z = 
the virial radius R vn - is defined as the radius of a 
sphere with overdensity of 360 of the average mat- 
ter density. The overdensity limit changes with 
redshift and asymptotically goes to 178 for high z. 
Different definitions are also found in the litera- 
ture. For example, the often used overdensity 200 
relative to the critical density gives mass Maooj 
which for Milky- Way-mass halos is about 1.2-1.3 
times smaller than M v ; r . The exact relation de- 
pends on halo concentration. 

Overall, there are about 10 million halos in Bol- 
shoi. Halo catalogs are complete for halos with 
V cilc > 50 km s- 1 (M vir w 1.5 x IO^/T^Mq). 
We track evolution of each halo in time using 
~180 stored snapshots. The time difference be- 
tween consecutive snapshots is rather small: ~(40- 
80) Myrs. 

4. Halo mass and concentration functions 

Throughout most of the paper we character- 
ize halos using their circular velocity. However, 
V^i rc tightly correlates with halo mass, as demon- 
strated in Figured! For distinct halos with M v - lr = 



10 12 - 1O 14 /i _1 M 90% of halos have their circu- 
lar velocities within 8% of the median value. Even 
99% of halos are within only 15-20%. The varia- 
tions are substantially larger for subhalos: 90% of 
subhalos with masses 10 11 — 10 13 ft- _1 M Q lie within 
20% of the median V C i Ic . On average, the circular 
velocity increases with mass. The V c i rc — M v - lr rela- 
tion depends on halo concentration c(Af v j r ), and, 
thus, studying this relation gives us a way to es- 
timate c(M v ; r ) without making fits to individual 
halo profiles. Any halo mass profile can be param- 
eterized as M(r) — Mof(r/ro), where Mo and ro 
are parameters with mass and radius units and the 
function f(x) is dimensionless. If x max is the di- 
mensionless radius corresponding to the maximum 
circular velocity, then we can write the following 
relation between th e maximum circula r velocity 
and the virial mass (jKlvpin et al.ll200ll ): 



f(x) 



G 



M, 



fip) ^max 
4-7T 



1/3 



U 3 

ln(l +x) - 

Rvir 



-A, 



M ■ 



R 



1/2 

M vir 1 / 3 ,(3) 

(4) 
(5) 

= 2.15. (6) 



Here A V i r is the overdensity limit that defines the 
virial radius; p cr and Om are the critical density 
and the contribution of matter to the average den- 
sity of the Universe. The first two equations are 
general relations for any density profile. Eqs. ([5]) 
and (j6|) are specific for NFW: r s is the characteris- 
tic radius of the NFW profile, which is the radius 
at which the logarithmic slope of the density is —2. 
At z = 0, for our cosmological model, A V j r = 360 
and f2 m = 0.27. Calculating the numerical factors 
in eqs.([3][6|) we get the following relation between 
virial mass, circular velocity and concentration at 
z = 0: 



v/ln(l + c)~ c/(l + c) • 

where mass M V - 1T is in units of H~ 1 Mq, circular 
velocity is in units of km s _1 . This relation gives 
us an opportunity to estimate halo concentration 
directly for given virial mass and circular velocity. 

Alternatively, one may skip the c(M) term 
and simply use power-law approximations for the 
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Fig. 4. — Dependence of maximum circular velocity V c - ac on halo mass for distinct halos at redshift z = (left 
panel) and redshift z = 3 (right panel). The circular velocity at any moment mostly scales as V c - nc oc M^ 1 ' 3 . 
The Figure shows deviations from this scaling law. The deviations are related with the halo concentration. 
Full curves on both plots show median V^rc and dashed curves show 90% limits. Dots represent individual 
halos with large masses. 



V c i rc — M v i r relation which give a good fit to nu- 
merical data. 
For distinct halos: 



Kirc(Afvir) = 2.8 x 10" 2 M V 



0.316 



and for subhalos: 



Kirc(M sub ) = 3.8 x 10" 2 M 



sub 



0.305 



(8) 



(9) 



Here velocities are in km s 1 and masses are in 
units of h~ 1 M ( 7 ) . 

Equations ((3] - E]) can be considered as equa- 
tions for halo concentration: for given z, M v - lr , and 
T4; rc one can solve them to find c. For quiet halos 
the result must be the same as what one gets from 
fitting halo density profiles with the NFW approx- 
imation: two independent parameters (in our case 
M v ; r and V^rc) uniquely define the density pro- 
file. However, by using only two parameters, we 
are prone to fluctuations. In order to reduce the 
effects of fluctuations we apply eqs. (j3] - [6]) to the 
median values of V^rc for each mass bin. For each 
mass bin with the average M V j r we find the me- 
dian circular velocity T4i rc (M v i r , z). We then solve 



equations ([3][6]) and get median halo concentra- 
tions c(Mvir, z). One can also find concentration 
for each halo and then take the median - result is 
the same because for a given mass the relation 
between c and V c i rc is monotonic. This proce- 
dure minimizes effects of fluctuations and gives 
the median halo concentrat ion for a given mass. 
Munoz-Cuartas et alJ (|2010h applied our method 
to their simulations and repr oduced result s of di- 
rect density profile fitting. iPrada et al.l (|201ll ) 



state that this method recovers results of halo con- 



centrations c(M) in MS-I (|Neto et al.1 120071 ) with 
deviations of less than 5 percent over the whole 
range of masses in the simulation. 

Figure [5] shows the concentrations for redshifts 
z = 0—5. For redshift zero we get the following 
approximation: 



c(M vir ) = 9.60 



M v 



-0.075 



for distinct halos. For subhalos we find: 

-0.12 

/ 1 / 1 \ 

c(M sub ) = 12 



sub 



10 12 h- l Mr. 



(10) 



(11) 



s 



10 - 




M vir (h-'M Q ) 



Fig. 5. — Evolution of concentration of distinct 
halos with redshift. The full curves and symbols 
show results of simulations. Analytical approxi- 
mations are shown as dashed curves. All the fits 
have the same functional form of eq. (fT"2")) with two 
free parameters. At low redshifts the halo concen- 
tration decreases with increasing mass. However, 
the trend changes at high redshifts when the con- 
centration is nearly flat and even has a tendency 
to slightly increase with mass. 

Subhalos are clearly more concentrated than 
distinct halos with the same mass, which is likely 
caused by tidal stripping. The differences between 
subhalos and distinct halos on average are not 
large: a 30% effect in halo concentration. 

We also study the evolution of distinct halo con- 
centration with redshift. Figure [5] shows the con- 
centrations for redshifts z = — 5. Results can be 
approximated using the following functions: 



c(M vir ,z) 



CoO) 



( M vi 



-0.075 



\M (z) 



o 

0.26' 



(12) 



where cq(z) and M (z) are two free factors for 
each z. Table [3] gives the parameters for this 
approximation at different redshifts. For conve- 
nience we also give the concentration for a virial 



Fig. 6. — Evolution of halo concentration for ha- 
los with two masses indicated on the plot. The 
dots show results of simulations. For the refer- 
ence the dashed lines show a power-law decline 
c oc (1 + -z) 1 - Concentrations do not change as 
fast as the law predicts. At low redshifts z < 2 
the decline in concentration is c oc 6 (dot-dashed 
curves), where S is the linear growth factor. At 
high redshifts the concentration flattens and then 
slightly increases with mass. For both masses the 



concentration reaches a minimum of c n 



4-4.5, 



but the minimum happens at different redshifts for 
different masses. The full curves are analytical fits 
with the functional form of eq. (JT3J) - 



mass 10 12 ft _1 M(3 and the minimum value of con- 
centration c m j n . The simulation box for the Bol- 
shoi is not large enough to find whether there is 
a minimum concentration for z < 0.5. For these 
epochs the Table gives the value of concentration 
at 10 15 /i _1 M Q as predicted by the analytical fits. 

The curves in Figure [5] look different for dif- 
ferent redshifts. Typically the concentration de- 
clines with redshift and the shape of the curves 
evolves. Another interesting result is that the 
high-z curves show that the halo concentration has 
an upturn: for the most massive halos the concen- 
tration increases with mass. In order to demon- 
strate this more clearly, we study in more detail 
the evolution with redshift of the halo concentra- 
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tion for halos with two masses: 3 x 10 11 M 
and 3 x 10 12 /i _1 Mq. Note that the masses are 
the same at different redshifts. So, this is not the 
evolution of the same halos. Figure [5] shows the 
results. Just as expected, in both cases at low 
redshifts the halo concentration declines with red- 
shift. The decline is not as steep as often assumed 
c oc (1 + it is significantly shallower even 

at low z. For z < 2 a power-law approximation 
c oc 8(z) is a much better fit, where 6(z) is the lin- 
ear growth factor. It is also a better approxima- 
tion because the evolution of concentration should 
be related with the growth of perturbations, not 
with the expansion of the universe. At larger z 
the concentration flattens and slightly increases at 
z > 3. The upturn is barely visible for the larger 
mass, but it is clearly seen for the 3 x 10 11 h" 1 Mq 
mass halos. These and other results show that 
the concentration in the upturn does not increase 
above c w 5 though it may be related with the 
finite box size of our simulation. There is also 
an indication that there is an absolute minimum 
of the concentration c min « 4 at high redshifts. 
Relaxed halot[f] show a slightly stronger upturn in- 
dicating that non-equilibrium effects are not the 
prime explanation for the increasing of the halo 
concentration. 

The following analytical approximations pro- 
vide fits for the evolution of concentrations for 
fixed masses as shown in Figure [6) 

c(M vir ,z) = c(M vir ,0) [5 4 / 3 {z) + K (S- 1 (z)-l)\ , 

(13) 

here 6(z) is the linear growth factor of fluctua- 
tions normalized to be 5(0) — 1 and k is a free 
parameter, which for the masses presented in the 
Figure is re = 0.084 for M = 3 x lO^h^M® and 
K = 0.135 for ten times more massive halos with 
M = 3 x 10 12 h- 1 M Q . 

It is interesting to compare these result s with 
other simulations. IZhao et all (|2003l l2009h were 
the first to find that the concentration flattens at 
large masses and at high redshifts. Their estimates 
of the minimum concentr ation are compatib le with 



Table 3: Parameters of fit eq.(fT2"j) for virial halo 
concentration 



our results. Figure 2 in IZhao et al.l (|2003l) shows 



Redshift 


CO 


Mo/h~ L M e 


Cmin 


c{W A h' L M Q ) 


0.0 


9.60 






9.60 


0.5 


7.08 


1.5 x 10 17 


5.2 


7.2 


1.0 


5.45 


2.5 x 10 15 


5.1 


5.8 


2.0 


3.67 


6.8 x 10 13 


4.6 


4.6 


3.0 


2.83 


6.3 x 10 12 


4.2 


4.4 


5.0 


2.34 


6.6 x 10 11 


4.0 


5.0 




10 10 10 11 10 12 10 13 10 14 10 15 
M . (h-!M n ) 

Fig. 7. — Mass function of distinct halos at dif- 
ferent redshifts (circles). Curves show the Sheth- 
Tormen approximation, which provides a very ac- 
curate fit at z = 0, but overpredicts the number 
of halos at higher redshifts. 

an upturn in concentration at z = 4. However, the 
results were noisy and inconclusive: the text does 
not even mention it. 

iMaccio et al. (2008) present results that can 
be directly compared with ours because they use 
the same definition of the virial radius and esti- 
mate masses within spherical regions. Their mod- 
els named WMAP5 have parameters that are very 
close to those of Bolshoi. T here is one potentia l 



3 Relaxed halos are defined as halos with offset parameter 
X Q ff < 0.07 and with spin parameter A < 0.1, where X g 
is the distance from the halo center to its center of mass in 
units of the virial radius. 



Maccio et ail (|2008h 



issue with their simulations, 
use a set of simulations with each simulation hav- 
ing a small number of particles and either a low 
resolution, if the box size is large, or very small 
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box, if the resolution is small. For all the halos 
in their simulations the approximation for concen- 
tration is c(Mvir) = 8.41(M vir /lO 12 /i- 1 M )- - 108 . 
Bolshoi definitely gives more concentrated halos. 
The largest difference is for cluster-size halos. For 
M vir = lO^h^Mp, our results give c = 5.7 



6 for halos with M v 



(1-10) x io"r% 



while iMaccio et all (|2008f ) predict c = 4.0 - a 
40% difference. The difference gets smaller for 
galaxy-size halos: 14% for M vir = 10 12 /i _1 Mg 
and 2% for M vir = 10 10 /i _1 Mq. Comparing re- 
sults for relaxed halos we find that the disa gree- 



give 
as com- 



ment is smaller. iMunoz-Cuartas et al. ( 2010h 
c(M vir ) = 9.8 for M vir = 1Q 12 h^M® 
pared with our results (also for relaxed halos) of 
c(M V i r ) = 10.1 - a 3% difference. For clusters with 
M vir = 10 15 h- 1 M Q the disagreement is 18%. 

We can also compare our results with those of 
MS-I (|Neto et al.l 120071) though MS-I has differ- 
ent cosmological parameters and a different power 
spectrum. Because the analysis of MS-I was done 
for the overdensity 200, we also made halo cat alogs 
for this definition of halos. iNeto et al.l (|2007l) give 
the following approximation for all halos: C200 = 
7.75(M 2 oo/10 12 /i- 1 Af Q )-° 11 . For halos in the 
Bolshoi simulation C200 = 7.2(M 200 /10 12 /i- 1 M©)" - 075 . 
Thus the MS-I C200 is 8% larger than the con- 



centrations in Bolshoi for M v 



10 



12 



with a small (~10%) difference for M v 
lO^/i^M©. For M vh 



h-'M e , 
10 14 - 
1 M . in Ihc M S- 



II and Aquarius simulations iBovlan-Kolchin et al 



(2009) give an even larger concentration of c V j r = 



12.9, which is 1.3 times larger than what we get 
from Bolshoi. Most of the differences are likely 
due to the larger amplitude of cosmological fluc- 
tuations in MS simulations because of the combi- 
nation of a larger erg and a steeper spectrum of 
fluctuations. 

The mass function of distinct halos is a clas- 



sical cosmological result (e.g., Warren et al. 2006: 


Reed et al. 2007: 


Tinker et al.l 2008: Maccio et al. 


20081: Reed et al. 


2009). Figure [7] presents the re- 



suits for the Bolshoi simulation together with the 
predictions of the Shet h-Tormen approximation 
(|Sheth fc Tormenl 120021 ST, see also Appendix 
B). We find that the ST approximation gives an 
accurate fit for z = with the deviations less 
than 10 percent for masses ranging from M v i T = 
5 x lO^i^M© to M vlr = 5 x 10 li h- 1 M Q . How- 
ever, the ST approximation overpredicts the abun- 
dance of halos at higher redshifts. For example, at 



the ST approximation gives a factor of 1.5 more 
halos as compared with the simulation. At red- 
shift ten the ST approximation gives a factor of 
ten more halos than what we find in the simula- 
tion. 

We introduce a simple correction factor which 
brings the analytical predictions much closer to 
the results of simulations. We find that the ST 
approximation multiplied by the following factor 
gives less than 10 percent deviations for masses 
5 x W 9 h- 1 M Q - 5 x 10 u h- 1 M Q and redshifts z = 
0-10: 



F(6) 



(5.501c5) 4 



(14) 



1 + (5.500(5) 4 ' 
where S is the linear growth rate factor normalized 
to unity at z — (see eq. (|B2jl ). 



Ou r results are in good agreement with lTinker et al 
(|2008l ). who present the evolution of the mass 
function for z = — 2.5 for halos defined using the 
spherical overdensity method. At redshift zero, 
their mass function for overdensity A = 200 rela- 
tive to the mean mass density is 20% above the ST 
approximation. This is expected because masses 
defined with spherical overdensity A = 200 are 
typically 10-20% larger than virial masses, which 
are us ed in our paper. Results in iTinker et al 



( 20081 ) together with our work indicate that at 
higher redshifts the mass function gets more and 
more below the ST approximation. In addition, 
the shape of the mass function in simulations gets 
steeper: there is a la r ger di sagreement at larger 



Tinker et al] (|2008| ) argue that this be- 



masses. 

havior indicates that the mass function is not "uni- 
versal" : it does not scale with the redshift only as 
a function of the amplitude of perturbations a(M) 
on scale M (see Appendix B for definitions). Our 
results extend this trend to redshifts at least as 
large as z = 10. Our res u lts als o qualitatively 
agree with ICohn fc White] (2008), who presem, 
spherical overdensity halo masses at z — 10. They 
also find substantially lower mass functions as 
compared with the ST approximation, though the 
differences with the approximation are somewhat 
smaller than what we find. 

These results are at odds with those ob- 
taine d with the Friends-Of-Friends (FOF) method 



Lukic et al.l 120071 iReed et al] 120071 120091: 



Cohn fc White! 120081 ). The FOF halo mass func 
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Fig. 8. — Velocity function of distinct halos at different redshifts. Left: Symbols and curves show the 
cumulative velocity function of the Bolshoi simulation. Error bars show Gaussian fluctuations. Right: We 
plot the product V 3 n(> V) of the cumulative velocity function and the circular velocity. Dotted curves show 
analytical approximations with the form given by V 3 n(> V) — Aexp(— [V/Vb] a ), which provide accurate 
fits to numerical results. 



tion scales very close to the "universal" a{M) 
behavior. The reason for the disagreement be- 
tween spherical overdensity and FOF halo find- 
ing methods is likely related with the fact that 
FOF has a tendency to link together structures 
before they become a part of a virialized halo. 
This happens more often with the rare most mas- 
sive halos, which have a tendency to be out of 
equilibrium and in the process of merging. As a 
result of this, FOF ma sses are artificially inflated. 
Cohn fc Whit! (2008) studied case-by-case some 
halos at z = 10 and conclude that in their simu- 
lations FOF assigned to halos almost twice more 
mass. Comparison of the ST predictions with the 
Bolshoi results shown in Figure [7] points to the 
difference of a factor of 2.5 in mass for z = 10. 
Because of the steep decline of the mass function, 
a factor of 2.5 increase in mass translates to a fac- 
tor of ten increase in the number-density of halos. 
This correction to the FOF masses must be taken 
into account when making any estimates of the 
frequency of appearance of high-z objects. 

In Appendix C we also directly compare FOF 
masses with those obtained with the BDM code. 



At z = 8.8 the FOF masses with the linking pa- 
rameter I — 0.20 were on average 1.4 times larger 
than the BDM masses. In addition, the spread of 
estimates was very large with FOF in many cases 
giving 3-5 times larger masses than BDM. Anal- 
ysis of individual cases shows that this happens 
because FOF links large fragments of filaments, 
not just an occasional neighboring halo. The sit- 
uation is different at z = 0. Here both BDM and 
FOF (7 = 0.17) give remarkably si milar results 



thoug h some spread is still present (j Tinker et al 



2008). We speculate that the difference in the be- 



havior at high and low z is related with the slope 
of the power spectrum of perturbations probed by 
halos at different redshifts. 

These differences between different definitions 
of masses and radii of halos indicate the inherent 
weakness of masses as halo properties: in the ab- 
sence of a well defined physical process responsible 
for halo formation, masses are defined somewhat 
arbitrarily. We know that halos do not form ac- 
cording to the often used top-hat model. We also 
know that the virial radius is ill-defined for non- 
isolated interacting objects. Nevertheless, we use 
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Fig. 9. — Parameters of the velocity function at 
different amplitudes of perturbations erg (z). Open 
circles show parameters found by fitting n(> V,z) 
at different redshifts. The curves are power-law 
fits given by eqs. (|15H17p . The top right panel 
shows the evolution of as with redshift as pre- 
dicted by the linear theory. Circles on the curve 
indicate the same moments as on the other three 
panels. 

one definition or another and we pay a price for 
this vagueness. These uncertainties in masses also 
motivate us to use another, - much better defined 
quantity - the maximum circular velocity. 

5. Halo velocity function 

The velocity function for distinct halos is shown 
in Figure [8] It declines very steeply with velocity. 
At small velocities the power slope is -3 with an 
exponential cutoff at large velocities. We find that 
at all redshifts the cumulative velocity function 
can be accurately approximated by the following 
expression: 



n(> V) = AV~ exp - 



V 



(15) 



where the parameters A, Vq, and a are functions 
of redshift. For z — we find 



.4 



1.82 x 10 4 (/i _1 Mpc/km s^ 1 ) 



Fig. 10. — The velocity function of satellites com- 
pared with the velocity function of distinct halos. 
The bottom panel shows the cumulative function 
of subhalos (bottom full curve) and distinct ha- 
los (top full curve). The circular velocity used for 
the plot is the peak over each halo's history. The 
dashed curves are analytical approximations. The 
top panel shows the ratio of the number of subha- 
los and distinct halos (full curve) and an analytical 
approximation for the ratio (dashed curve). 



a = 2.5, 

Vo = 800 km s- 1 . 



(16) 



The evolution of the parameters should not di- 
rectly depend on the redshift, but on the ampli- 
tude of perturbations. Indeed, when we plot the 
parameters as functions of as(z) as predicted by 
the linear theory at different redshifts, the func- 
tions are very close to power-laws as demonstrated 
by Figure [51 We find the following fits to the pa- 
rameters: 



.4 
a 

Vo 



1.52 x 10 4 ct 8 3/4 (z)(/i" 1 Mpc/km s" 1 ) 



-l\-3 



8 

4/3/ 



1 + 2.15cts /j (z), 



(17) 



3300 



1 + 2.5ct|(z) 



-km s 



6. Abundance of Subhalos 

Figure 1101 shows the cumulative velocity func- 
tion n(> V c i rc ) of all subhalos in Bolshoi regardless 
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of the circular velocity of their host halos. We use 
maximum circular velocities over the whole evolu- 
tion for both subhalos and distinct halos. The top 
panel shows the ratio of the number of subhalos to 
the number of distinct halos with the same limit 
on the circular velocity. Note that for given V c i IC 
most of the halos are distinct. This may sound a 
bit counter intuitive. Because each distinct halo 
has many subhalos, one would naively expect that 
there are many more satellites as compared with 
distinct halos. This is not true. The number of 
satellites is large, but most are small. When we 
count small distinct halos, their number increases 
very fast and we always end up with more distinct 
halos at a given circular velocity. The abundance 
of subhalos can be approximated with the same 
function given by cq. (p~5|) as for the distinct halos. 
However, the parameters of the approximation are 
different. For subhalos that exist at z — and for 
which we use peak velocities over their history of 
evolution, we find: 

A = 6.2 x 10 3 (^ 1 Af Q /km s" 1 )" 3 

a = 2.2, V = 48010ns- 1 . (18) 

The remarkable similarity of the shapes of the ve- 
locity functions of halos and subhalos suggests a 
simple interpretation for the difference in their pa- 
rameters: subhalos were typically accreted at the 
epoch when the velocity function of distinct halos 
had the same parameters a and Vq as the velocity 
function of subhalos at present. For the parame- 
ters given by eqs. ()17|18[) we get a typical accretion 
redshift z acc sa 1. 

In order to study statistics of subhalos belong- 
ing to different parent halos we split our sample 
of distinct halos into subsamples with different 
ranges of circular velocities. For each subhalo we 
use either its z = circular velocity or the peak 
value over its entire evolution. Figure [11] shows 
both the present day and the peak velocity distri- 
bution functions for distinct halos with masses and 
velocities ranging from galaxy-size halos to clus- 
ters of galaxies. The average circular velocities for 
each bin presented in the figure are (from bottom 
to top): Vhost = (163,190,235,340,470,677,936) 
km s . The number of halos in each bin varies 
from 200 for the most massive halos (Vh os t = 
800 - 1200 km s" 1 ) to 30000 for the least massive 
halos with V host = 160 - 180 km s -1 . The in- 




Fig. 11. — The cumulative velocity function of 
satellites for host halos with different maximum 



to 



circular velocities ranging from « 150 km s" 
« 1000 km s _1 from bottom to top. The bottom 
panel uses the velocities of subhalos at redshift 
z = 0. The top panel uses peak circular velocities 
over the history of each subhalo. The dashed lines 
show power-laws with the slope —3. The abun- 
dance of subhalos increases with increasing host 
halo mass. 



crease in the abundance of substructure for more 
hosts is consistent with the results of 



massive 



Gao et al.1 (|2004l ). who give a factor of 2.0-2.5 in- 
crease for host halos from mass ~ 2.5xl0 12 h~ 1 M^ l 
to ~ 10 15 h' 1 Mr*. Ivan den Bosch et all ([2005h : 
Taylor fc Babull (|2005l) : IZentner et all(|2005l ) came 
to similar conclusions using their (semi) analytic 
models. 

In order to more accurately measure the de- 
pendence of the abundance of subhalos on the 
circular velocity of the host halo, we analyze the 
number of satellites with circular velocities larger 
than 0.3 of the circular velocity of their hosts: 
V sa _t > 0.3 Vhost- This approximately corresponds 
to the mass ratio of M sat /M host ~ 0.3 3 « 0.027. 
The threshold of 0.3 is a compromise between the 
statistics of satellites and the numerical resolu- 
tion of the simulation. Figure Q21 shows the num- 
ber of satellites iV .3(Vhost) for hosts ranging from 
Vhost ~150 km s" 1 to -1000 km s _1 . The num- 
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Fig. 12. — Dependence of the number of subhalos 
on the circular velocity of their hosts. Here we 
count all subhalos with circular velocities larger 
than 0.3 of their host circular velocity. The bot- 
tom panel shows results for y circ estimated at 
z = and the top panel is for the peak V C i IC over 
the history of each subhalo. For z = circular 

1 /2 

velocities the abundance scales as ^ hos t (dashed 
curve). For peak circular velocities the number 
of subhalos is larger and the scaling is steeper: 
N oc V^l*t ( lun curve). For comparison, the 
dashed curve is the same as on the bottom panel. 



ber of satellites scales as a power-law A*o.3 oc VCL t 
with the slope 7 depending on how the circular 
velocity is estimated. For the z = velocities the 
slope is 7 = 1/2, and it is larger for the peak ve- 
locities: 7 = 2/3. 

There is an indication that the dependence of 
the cumulative number of satellites on their circu- 
lar velocity N{> V s&t ) gets slightly shallower for 
more massive host halos. Figure [T31 illustrates the 
point. Here we study the most massive (but also 
rare) halos. Again, the number of satellites is ap- 
proximated by a power-law. However, the slope 
is about —2.75, which is somewhat shallower than 
the slope —3 found for smaller host halos. 

We compare some of our results with the 



Via L actea-II (VL-II) simulation (jDiemand et al 
2008). We do not use the published results be- 



V at (km/s) 



Fig. 13. — Velocity function of subhalos for 
40 most massive clusters with average M vir = 
3.2 x lQ li h~ 1 M Q and circular velocity V C i TC = 
1100 km s -1 . The velocity function is nearly a 
power law with the slope -2.75. 



cause the analysis of VL-II was done using over- 
density 180 relative to the matter density, which 
gives a larger radius for halos as compared with 
the virial radius. We use the halo catalog of 
VL-II, which lists coordinates and circular ve- 
locities of individual halos. We also use pub- 
lished parameterization of the dark matter den- 
sity in order to estimate the virial radius of VL- 
II. When comparing with VL-II, we select ha- 
los in Bolshoi in a narrow range of circular ve- 
locities V^i rc = 195 — 205 km s _1 . There are 
4960 of those with the average virial mass of 
M vir = 1.26 x lO 12 /& _x Af , which is close to 
the virial mass M vir = 1.3 x 1O 12 /i _1 M of Via 
Lactea II. Figure Q3] presents results of the ve- 
locity function of subhalos in those host halos. 
The dashed line in the figure is the power law 
N(> V) = (V/61km s -1 ) -3 , which gives a good 
fit to the data for a wide range of circular ve- 
locities from 4 km s _1 to 100 km s — 1 . Bolshoi 
has slightly more subhalos by about 10%. This 
is a small difference and it goes in line with ex- 
pectation that a smaller normalization of cosmo- 
logical fluctuations gives fewer subhalos. In the 
same vein, the Aquarius simulations have an even 



15 




Fig. 14. — Comparison of satellite velocity func- 
tions in the Via Lactea-II and Bolshoi simula- 
tions for host halos with V C i rc = 200 kms/s and 
M vir « 1.3 x I0 12 h- 1 M Q . The dashed line is a 
power law with slope —3 which provides an excel- 
lent fit to both simulations. In both simulations 
satellites are found inside a sphere with virial ra- 
dius i? V i r . 

higher (by 30% as compared to VL-II) number 
of subhalos, probably because of an even larger 
amplitude. 

Summarizing all the results, we conclude that 
the cumulative velocity function of z — subhalos 
can be reasonably accurately approximated by the 
power law: 



v host x ! 



N(>x) = 1.7x10" 

x = V suh /V host , x < 0.7, 



(19) 
(20) 



where the circular velocity of the host is given 
in units of km s . Again, these results arc 
broadly consisten t with the iV-body simulations of 
Gao et al.l (120041) and with the s e mianalytic mod- 
els of van den Bosch etail (l2005l ); lTavlor fc Babul 
(|2005l) and lZentner et alT(|2005[ ). For peak circu- 
lar velocities we obtain: 



N(> x) = 9.0 x lQ~Xo^ 



(21) 



Fig. 15. — Density profiles for galaxy-size halos 
with V C i IC ~ 200 km s -1 . The full curve and cir- 
cles are the dark matter density and the number- 
density of satellites with V CU c > 4 km s _1 in 
the Via Lactea-II simulation normalized to the 
average (number-)density for each component re- 
spectively. The satellites have nearly the same 
overdensity as the dark matter for radii R = 
(0.3 — 2)i? v ; r . The number-density of satellites 
falls below the dark matter at smaller radii. The 
dashed curve is the number-density of satellites 
with V^irc > 80 km s _1 found at z — in the 
Bolshoi simulation for host halos selected to have 
the same circular velocity as Via Lactea-II. In the 
outer regions with R — (0.5— 1.5)i? v i r the satellites 
follow the dark matter very closely. In the inner 
regions the Bolshoi results are 20-30% below the 
much higher resolution simulation Via Lactea-II, 
presumably because of numerical effects. 



7. The spatial distribution of satellites 

The spatial distribution of satellites has nu- 
merous astrophysical applications. Among oth- 
ers, these include the survivability of dark mat- 
ter subhalos (e.g.jMoore et alJll996t iKlypin et al 



19991 : IColm et al.l [1999J), the pote ntial annihila- 
tion signal of dark matter (e.g., K uhlen et al 
20081 : ISpringel et al.l l2008t lAndoll2009h . and mo- 
tion of satellites as a probe for masses of iso- 
lated galaxies and groups (e.g., Prada et al. 20031: 
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Klypin fc Pradal 120091: iMore et all l2009h . The 
relative abundance of satellites and dark mat- 
ter is a form of bias. Thus, studying the distri- 
bution of satellites in simulations sheds light on 
the physics of bias and, thus, on the formation 
of dwarf galaxies. There is an additional reason 
to study the satellites in the Bolshoi simulation: 
comparison with high resolution simulations such 
as Via Lactea gives an additional test on resolu- 
tion effects and provides limits of the applicability 
of the simulation. 

The spatial distribution of satell ites has been 
extensively studied in simulations (IGhigna et al 



extensively studied m simulations (umgna et al. 
1998[l2000l:lNagai fc Kravtsoyl2005[lDiemand et"al1 
200a ISpringel et all l2008t lAngulo et all 120091 
One of the main issues regarding satellites is to 
what degree their distribution is more extended 
than that of the dark matter. As a small halo 
falls into the gravitational potential of a larger 
halo, it experiences tidal stripping and dynamical 
friction. It may also experience interaction with 
other subhalos before and during mfall. Tidal 
stripping reduces the mass of subhalos, resulting 
in a very strong radial bias: subhalos selected by 
mass have relatively low number-density in the 
central region of their hosts. However, stripping 
affects much less the central parts of subhalos. 
This is why the distribution of subhalos is more 
concen trated when selected by their circular ve- 
locity (jNagai fc Kravtsovl 120051 ) . Depending on 
the mass and concentration of the subhalo and on 
its trajectory, the role of different physical effects 
may vary. Interplay of these processes results in 
a complicated picture of the distribution of the 
satellites. Numerical effects add to the complex- 
ity of the situation: it is a challenge to preserve 
and to identify small subhalos through the whole 
history of evolution of the Universe. 

The traditional way of displaying results is to 
normalize both the dark matter and satellites to 
the average mass inside the virial radius. When 
presented in this way, results routinely show that 
there are relatively more s atellites in the outer 
parts of halos. For example, Angulo et al. (2009) 
find that independently of subhalo mass, subha- 
los are a factor of two more abundant than dark 
matter around the virial radius of their hosts. If 
that were correct, this would imply some kind 
of physical mechanism to produce more satellites 
outside of the virial radius, a far-reaching conclu- 




0.4 0.6 0.8 1 
R/R . 

Fig. 16. — Comparison of satellites (dashed 
curves) and dark matter density (full curves) ra- 
dial distributions for halos with M V i r = 5 x 
10 12 h~ 1 M Q (bottom panel) and A/ vir = 3 x 
IO^-^q (top panel). Halos were selected to 
be isolated: no other larger halo within radius 
2R vir . Subhalos fallow very closely the dark mat- 
ter density in the outer regions R > 0.5i? vir . In 
the central regions of the halos the overdensity of 
satellites is smaller than that of dark matter. 



sion. However, below we show that this not cor- 
rect. 

The main issue here is the normalization. If 
a simulation includes only a small region, a typi- 
cal setup for modern high-resolution simulations, 
there is no sensible way to normalize satellites. 
This is not the case given the statistics of Bol- 
shoi. We can reliably normalize the abundance of 
even halos with small mass and circular velocities. 
The velocity function of all distinct halos is given 
by eq. (jTB"j) . There are 20-25% subhalos at given 
circular velocity. Using the velocity function and 
the fraction of subhalos from Bolshoi, we estimate 
the average number of all halos with any given 
Vcitc. These estimates are used to normalize the 
number-density profile of subhalos in VL II pre- 
sented in Figure [15] A comparison with the dark 
matter profile is quite interesting: there is very lit- 
tle bias in the Via Lactea distribution of subhalos 
for radii R — (0.3 — 2)i? v ; r . Subhalos are not more 
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extended as compared with the dark matter. How- 
ever, Via Lactea-II is just one halo and there may 
be some effects related with cosmic variance. We 
use halos in Bolshoi that have similar characteris- 
tics as Via Lactea-II: very isolated halos (no equal 
mass halo inside 3i? v ir) and circular velocities in 
the range V c i IC = (200 — 220) km s _1 with corre- 
sponding masses M vir = (1.3— 1.5) x 10 12 h~ 1 M Q . 
For Bolshoi halos we find a small ~10% antibias. 
In the inner regions (R < 0.5i? v ir) the number 
density of satellites goes below the results of Via 
Lactea, but the difference is not large: 20-30%. 
Some of the differences with Via Lactea-II may be 
real because subhalos in Bolshoi are more massive, 
and, thus, they must experience stronger dynami- 
cal friction. However, it is more likely that most of 
the differences are numerical: after all, Bolshoi has 
substantially worse resolution than Via Lactea-II. 
Regardless of the cause of those small differences, 
it is quite remarkable that simulations with five 
orders of magnitude difference in mass resolution 
produce results that deviate only by 10-20%. 

Comparison with Via Lactea-II is difficult be- 
cause for these masses (M v i r s» 10 12 /j _1 Mq) Bol- 
shoi has only a few subhalos per each host. In or- 
der to have a better picture of the spatial distribu- 
tion of satellites, we study more massive halos for 
which our resolution is relatively better. Again, 
we select isolated halos: those with no larger halo 
within twice the virial radius. Figure [TBI shows the 
results for hosts with very different masses. The 
top panel shows results for 82 halos with circu- 
lar velocities in the range 900-1100 km s _1 (av- 
erage virial mass 2.5 x 10 14 /i _1 Mq). The bot- 
tom panel shows 2200 halos with V^ rc = 280 — 
300) km s" 1 and average M vir = 5x 10 12 /i _1 M Q . 
Results are very similar for such different host ha- 
los: satellites follow the dark matter very closely 
for radii R — (0.5 — 2)i? V i r with possible small 
(~10%) antibias. In the central region the sub- 
halo abundance goes below the dark matter by a 
factor of 2-2.5 at R = 0.2i? vir . 

It is interesting to compare the Bolshoi results 
with iNagai fc Kravtsovl (|2005h . who present pro- 
files for 8 clusters with almost the same masses 
as in the top panel of our Figure 1161 If we 
change their normalization for subhalos selected 
by present-day circular velocity (their Figure 3) in 
such a way that the dark matter profile matches 
the overdensity of satellites at the virial radius (we 



use a factor of 0.8 to match our definition of virial 
radius), then there is excellent agreement with 
Bolshoi, with both simulations giving the ratio of 
the dark matter density to the number-density of 
satellites ~ 2 at R = 0.2i? V i r . 

8. Conclusions 

Using the large halo statistics and high reso- 
lution of the Bolshoi simulation we study numer- 
ous properties of halos and subhalos. We present 
accurate analytical approximations for such char- 
acteristics as the halo and subhalo abundances 
and concentrations, the velocity functions, and the 
number-density profiles of subhalos. Detailed dis- 
cussions of different statistics have already been 
given in relevant sections of the text. Here we 
present a short summary of our main conclusions. 

Velocity function. Our main property of ha- 
los is their maximum circular velocity V c i TC . As 
compared with virial masses, circular velocities 
are better quantities to characterize the physical 
parameters of the central regions of dark matter 
halos. As such, they are better quantities to re- 
late the dark matter halos and galaxies, which 
they host. We present the halo velocity functions 
at different redshifts and show that they can be 
accurately described by eqs. (115H17[) . The halo 
circular velocity function n(> V) declines as a 
power-law V~ 3 at small velocities and has a quasi- 
exponential cutoff at large circular velocities. 

Mass function of disti nct halos. We find that 



the ST approximation (jSheth fc Tormenl 120021 ) 
gives an accurate fit to the redshift-zero mass func- 
tion: errors are less than 10% for masses in the 
range 5 x lO lo /i _1 M - 5 x lO 14 fr _1 M . However, 
the approximation ovcrpredicts the halo abun- 
dance at higher redshifts and gives a factor of ten 
more halos than the Bolshoi simulation at z = 10. 
The correction factor eq. (fT4|) brings the accuracy 
of the approximation back to the ~ 10% level for 
redshifts z = — 10. It also breaks the universal- 
ity of the fit: the mass function cannot be written 
as a function of only the rms fluctuation a(M) 
on mass scales M. These results depend on how 
halos are defined with the Friends-Of-Friends algo- 
rithm giving different answers than the spherical 
overdensity method. See Appendix C for details. 

Concentrations of halos. The halo concentra- 
tion c(M y ij,z) appears to be more complex than 
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previously envisioned. For a given redshift z, the 
concentration first declines with increasing mass. 
Then it flattens-out and reaches a minimum of 
Cmin ~ 4 — 5 with the value of the minimum chang- 
ing with redshift. At even larger masses c(M v j r ) 
starts to slightly increase. This "up-turn" in the 
concentration is a weak feature: the change in 
concentration is only 20%. Moreover, it cannot 
be detected at low redshifts, z < 0.5. If our es- 
timates are correct, at z = the upturn should 
start at masses about M V - 1T ~ lO 18 ^ 1 M & - clus- 
ters this massive do not exist. However, at z > 2 
the upturn is visible at the very massive tail of 
the mass function. It is not clear what causes 
the upturn. The upturn is even stronger for re- 
laxed halos, which indicates that non-equilibrium 
effects cannot be the reason for the upturn. At 
large redshifts the halos that show the upturn are 
very rare: their mass is much larger than the char- 
acteristic mass of halos existing at that time. 
Most of them likely experience a very fast growth. 
They also represent very high-cr peaks of the den- 
sity field. It is known that the statistics of rare 
peaks are different from th ose of more "normal" 
peaks ( Bardeen et al. 19861 ). One may speculate 
that this may result in a change in halo concen- 
tration. More ext ensive analys i s of h alo concen- 
trations is given in lPrada et al. ( 2011 ). 

Subhalo abundance. The cumulative abundance 
of satellites is a power-law with a steep slope: 
N(> V) oc V~ 3 . Combing our re sults with those 



of th e Via Lactea-II simulation (|Diemand et al 
2008), we show that the power-law extends at 
least from 4 km s^ 1 to 150 km s _1 for Milky- Way- 
mass halos and yields the correct abundance of 
large satellites such as the Large Magellanic Cloud 
( Busha et al. 2010h . The abundance of satellites 
increases with the circular velocity and mass of 

1 /2 

the host halo as N oc V h ^ st . For example, this 
means that in relative units a cluster of galaxies 
has 2.5 times more satellites than the Milky Milky 
Way, in good agreement with previous numerical 



results (|Gao et alJl2004h . Eqs. (f20ll2T]) give ap- 



proximations for the abundance of satellites. 

Subhalo number- density distribution. One of 
the main issues here is the number-density of satel- 
lites relative to the dark matter in the outer re- 
gions of halos. Some previous simulations indi- 
cated substantial (a factor of two) overabundance 
of satellites around the virial radius. We do not 



confirm this conclusion. Our re-analysis of the Via 
Lactea-II simulation as well as the results from 
the Bolshoi simulation unambiguously show that 
there is no overabundance of satellites. In the 
Via Lactea-II simulation the satellites follow very 
closely the distribution of dark matter for radii 
R = (0.3 — 2)i? V j r . In the Bolshoi simulation we 
find a small (10%) antibias at the virial radius. 
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A. APPENDIX: Bound Density Maxima halofinder 

We u se a parallel (MPI+OpenM P) version of the Bound-Density-Maxima algorith m to identify halos i n 
Bolshoi ( Klypin fc Holtzman 19971 ). For detailed comparison with other halofinders see lKnebe et al. ( 2011 ). 



The code detects both distinct halos and subhalos. The code locates maxima of density in the distribution 
of particles, removes unbound particles, and provides several statistics for halos including virial mass and 
radius, as well as maximum circular velocity. The parameters of the BDM halo finder were set such that 
the density maxima are not allowed to be closer than 10/i _1 kpc. We keep only the more massive density 
maximum if that happens. This is mostly done to save computer time. It is also consistent with the 
force resolution of Bolshoi. Halo catalogs obtained with a smaller minimum separation of 7.5 /i _1 kpc did not 
include more halos. 

Removal of unbound particles is done iteratively. It goes in steps: 

1. Find the bulk velocity of a halo: the velocity with which the halo moves in space. The rms velocities 
of individual particles are later found relative to this velocity. We use the central region of the halo 
(the 30 particles closest to the halo center) to find the bulk velocity. 

2. Find the halo radius: the minimum of the virial radius and the radius of the declining part of the 
density profile (radius of the density minimum, if it exists). 

3. Find the rms velocity of dark matter particles and the circular velocity profile. Estimate the halo 
concentration. 

4. Find the escape velocity as a function of radius and remove particles that exceed the escape velocity. 
Use only bound particles for the next iteration. 

The whole procedure (steps 1-4) is repeated 4 times. If the mass or radius of a halo are too small (too few 
particles), the density maximum is removed from the list of halo candidates. 

If two halos (a) are separated by less than one virial radius, (b) have masses that differ by less than a 
factor of 1.5, and (c) have a relative velocity less than 0.15 of the rms velocity of dark matter particles 
inside the halos, we remove the smaller halo and keep only the larger one. This is done to remove a defect 
of halo-finding where the same halo is identified more than once. This removal of "duplicates" (halos with 
nearly the same mass, position, and velocity) happens only during major merger events when instead of two 
merging nearly equal-mass halos the halo finder sometimes finds 3-5 halos. Unfortunately, this also has the 
side effect of removing one of the major merger halos. This is a relatively rare event and it affects only the 
very tip of the subhalo velocity function. 

We use the virial mass definition M v ; r that follows from the top-hat model in an expanding Universe 
with a cosmological constant. We define the virial radius i? v ir of halos as the radius within which the mean 
density is the virial overdensity times the mean universal matter density p m = f^MPcrit at that redshift. 
Thus, the virial mass is given by 

M vir = —A virPm R 3 vir . (Al) 

Eq. (IBip gives an analytical approximation for A vir . For our set of cosmological parameters, at z — the 
virial radius R v - lT is defined as the radius of a sphere with overdensity of 360 times the average matter 
density. The overdensity limit changes with redshift and asymptotically approaches 178 for high z. 

Overall, there are about 10 million halos in Bolshoi (8.8 M at z = 0, 12.3 M at z = 2, 4.8 M at z = 5). Halo 
catalogs are complete for halos with U C i rc > 50 km s _1 (M vir w 1.5 x lO 10 ^ -1 M ). We do post-processing of 
identified halos. In particular, for distinct halos we find their properties (e.g., mass, circular velocity, density 
profiles) without removing unbound particles. For most, but not all, halos it makes little difference. For 
example, the differences in circular velocities are less than a percent for halos with and without unbound 
particles. Differences in mass can be a few percent depending on halo mass and on environment. 



4 We keep the peak that has the largest mass inside a sphere of radius 10 h 1 kpc 



22 




1 3 3 4567B9 10 

1+z 

Fig. 17. — Fraction of z = halos tracked to a given redshift for halos with different circular velocities at 
redshift zero. 

In order to track the evolution of halos over time, we find and store the 50 most bound particles (fewer, 
if the halo does not have 50 particles). Together with other parameters of the halo (coordinates, velocities, 
virial mass and circular velocity) the information on most bound particles is used to identify the same halos 
at different moments of time. The procedure of halo tracking starts at z = and goes back in time. The 
final result is the history (track) of the major progenitor of a given halo. The halo track may be lost at some 
high redshift when the halo either becomes too small to be detected or the tracking algorithm fails to find it. 
A new halo track may be initiated at some redshift if there is a halo for which there was no track at previous 
snapshots (smaller redshifts). This happens when a halo merges and gets absorbed by another halo. 

With ~180 snapshots stored, the time difference between consecutive snapshots is rather small. For 
example, the snapshot before the z = snapshot has z — 0.0027 with a time difference of 37 Myrs. The 
difference in time between snapshots stays on nearly the same level (42-46 Myrs) until z = 0.23 when it 
becomes twice as large. We start with z — halos and identify them in the previous snapshot. If a halo is 
not found at that snapshot, we try the next one. Altogether, we may try 6 snapshots. Typically, 95% of 
halos are found in the previous snapshot, an additional 2-3% in the next one and ~1% in even earlier ones. 
Overall, about (0.2-0.3)% of halos cannot be tracked at any given snapshot: they are either lost because 
their progenitor gets too small or because of numerical problems. The number depends on the redshift and 
on halo mass. Figure [T7] shows the fraction of halos tracked to given redshift for halos that exist at z = 0. 
More massive halos are tracked to larger redshifts. Half of all halos with I4i rc = 50 km s _1 are tracked to 
z = 4 and half of all halos with V c i TC — 200 km s^ 1 are tracked to z = 7. 

B. APPENDIX: Auxiliary approximations 

For completeness, here we present some approximations used in the text. For the family of fiat cosmologies 
(£Im + = 1) an accurate approximation for the value of the virial overdensity A v ; r is given by the analytic 
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formula (jBrvan fc NormarJ ll998): 

A vir = (18?r 2 + 82a; - 39x 2 )/n(z), 

where Q(z) = p m (z) / Pent an d x = Q(z) — 1. 

The linear growth-rate function 8(a) used in ag(a) is defined as 

8(a) = D(a)/D(l), 

where a = 1/(1 + z) is the expansion parameter and D(a) is: 



D(a) = - ( ^ M, ° 
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where f?M o an d ^a,o are the density contributions of matter and the cosmological constant at z = 
respectively. For 17m > 0-1 the growth rate facto r D(a) can be accurately approximated by the following 
expressions (jLahav et a,l.lll99lb ICarroll et alJll992h : 
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ft\(a) = 1 — Qm(«)- 

For fiM.o = 0.27 the error of this approximation is less than 7 x 10~ 4 . 

The Sheth-Tormen approximation (jSheth fc Tormenl 120021) for the distinct halos mass function can be 
written in the following form: 

dn da(M) 
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where M is the halo virial mass and 
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A = 0.322, 6 = 0.707. 
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Here P(k) is the power spectrum of perturbations, and W(k, M) is the Fourier transform of the real-space 
top-hat filter corresponding to a sphere of mass M. For the cosmological parameters of the Bolshoi simulation 
the rms density fluctuation a(M) can be approximated by the following expression: 

16.9y 041 
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The accuracy of this approximation is better than 2 percent for masses M > l0 7 h~ 1 M@. 
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Fig. 18. — Left: The halo mass function at redshift z — 8.8. The full curve shows the Sheth-Tormen (ST) 
approximation. Open circles show FOF halos identified using the linking length I = 0.20. The spherical 
overdensity halos are represented by solid circles. The ST approximation overpredicts the FOF (SO) masses 
by a factor 1.15 (1.9). Right: The distribution of mass around 9 massive halos (Mfof ~ 10 11 ft, -1 Mq) at 
redshift z — 8.8. Each panel shows 1/2 of the dark matter particles in cubes of 1 hr 1 Mpc size. The center 
of each cube is the exact position of the center of mass of the corresponding FOF halo. The effective radius 
of each FOF halo in the plots is 150 — 200ft, _1 kpc. Circles indicate distinct halos and subhalos identified 
by the spherical overdensity algorithm BDM. The radius of each circle is equal to the virial radius of the 
halo. The numbers in the top-left corner of each panel show the ratio of FOF mass to that of SO. Panels 
(a, c, g) show relatively good cases when the center of a halo in the simulation is close to the center of a 
FOF-detected halo. Panel (e) shows a major-merger: FOF linked the two halos together. In panels (b, d, f, 
h, i) FOF linked together halos which formed long and dense filaments. 



C. APPENDIX: FOF and SO masses 

In order to clarify the situation with the difference between the results on the halo mass function in the 
Bolshoi simulation and in the ST approximation at high redshifts, we make more detailed analysis of halos 
at redshift z = 8.8. We also study results obtained using the FOF halofinder with three linking-lengths: 
I = 0.17,0.20, and 0.23. We start by considering only the most massive halos with masses larger than 
10 11 /i _1 M Q . Each of those halos should have more than 700 particles. The spherical overdensity algorithm 
(BDM) identified 55 halos above this mass threshold. The FOF found 121, 255, and 602 halos with linking 
lengths I — 0.17, 0.20, and 0.23 above the same mass threshold. It is clear that FOF gives significantly higher 
masses. It is also very sensitive to the particular choice of the linking length. 

At z — 8.8 the ST approximation predicts a factor of 4 - 6 more halos as compared with what we find 
using the spherical overdensity algorithm. Because FOF with I = 0.20 gives about five time s more halos, it 
makes a v ery good match to the ST approximation. This is consistent with the results of Cohn fc White! 



(2008) and Reed et ail ( 20091 ). The left panel in Figure [TBI shows the mass functions at z = 8.8. At all masses 



FOF with / = 0.20 is well above the SO results and is close to the ST predictions. 

However, FOF results are very misleading. We compare the SO halos (as found by the BDM code) with 
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Fig. 19. — Left panels: Ratio of masses for the same halos identified by FOF with I = 0.20 and by the 
spherical overdensity BDM halofmders at redshift z = 8.8. On average, the FOF halos have masses 1.4 times 
larger than those obtained using SO. In addition, there is a significant spread in the mass ratios. Right 
panels: The same for z = with a linking-lcngth I = 0.17. Top panels show the distribution of mass ratios 
Mso/-^fof with halos selected by the SO mass. Bottom panels show the mass ratios for individual halos. 

the FOF halos found using I = 0.20. For halos with more than 100 particles both algorithms find essentially 
the same distinct halos, but FOF typically assigns larger masses to the same halos. The right panel in 
Figure [18] illustrates the point. There is a large variety of situations. We typically find that when there is 
a well-defined halo center and the halo dominates its environment, both the FOF and the SO masses are 
reasonably consistent (e.g., panel (a)). However, FOF has a tendency to link fragments of long filaments. In 
such cases the formal center of the FOF halo may not even be found in a large halo. Surprisingly, there are 
many of those long filaments at the high redshifts. 

Figure [19] presents statistics for the ratios of FOF and SO masses. Left panels show the most massive 
17000 halos with SO masses larger than 10 10 /i _1 M Q at z — 8.8. There is a large spread of masses and on 
average FOF masses are 1.4 times larger than the SO counterparts. We made the same analysis for the 
most massive 10000 halos with the SO masses larger than 5 x 10 12 ft. -1 M Q at z = using I = 0.17 for FOF. 
Remarkably, both halofmders produce similar results - a big contrast with high redshifts. Overall, there is 
a small offset in the mass ratios with SO producing 1.05 times larger masses. However, the difference is 
remarkably small. 

Thus, as we stated in §4, FOF halo masses are similar to SO ones at low redshifts, but systematically 
larger at high redshifts. 
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